/** \file
   \brief Vectorized functions.

   This file implements 
 **/


#ifndef _MMSL_VEC__HPP_
#define _MMSL_VEC__HPP_

namespace mmsl{




  template<class T, class itype=int>
  inline void 
  xeqa(const itype n, const T a, T *x, const itype incx)
  {
    itype i;
  #pragma omp parallel for default(none)		\
    shared(x) private(i)
    for (i = 0; i < n; ++i)
      x[i*incx] = a;
  }
  
  
  template<class T, class itype=int>
  inline void 
  xeqa(const itype n, const T a, T *x)
  {
    itype i;
#pragma omp parallel for default(none)		\
  shared(x) private(i)
    for (i = 0; i < n; ++i)
      x[i] = a;
  }
  
  
  template<class T, class itype=int>
  inline void 
  yeqapbx(const itype n, const T a, const T b, const T *x, const itype incx, T *y, const itype incy)
  {
    itype i;
   #pragma omp parallel for default(none)		\
     shared(x,y) private(i)
    for (i = 0; i < n; ++i)
      y[incy*i] = a + b * x[incx*i];
    
  }
  
  
  template<class T, class itype=int>
  inline void 
  yeqapbx(const itype n, const T a, const T b, const T *x, T *y)
  {
    itype i;
   #pragma omp parallel for default(none)		\
     shared(x,y) private(i)
    for (i = 0; i < n; ++i)
      y[i] = a + b*x[i];
    
  }
  
  
  template<class T, class itype=int>
  inline void
  yeqx(const itype n, const T *x, const itype incx, T *y, const itype incy)
  {
    itype i;
   #pragma omp parallel for default(none)		\
     shared(x,y) private(i)
    for (i = 0; i < n; ++i)
      y[incy*i] = x[incx*i];
    
  }
  
  template<class T, class itype=int>
  inline void 
  yeqx(const itype n, const T *x, T *y)
  {
    itype i;
   #pragma omp parallel for default(none)		\
     shared(x,y) private(i)
    for (i = 0; i < n; ++i)
      y[i] = x[i];
  }
  
  template<class T, class itype=int>
  inline void 
  yeqax(const itype n, const T a, const T *x, const itype incx, T *y, const itype incy)
  {
    itype i;
   #pragma omp parallel for default(none)		\
     shared(x,y) private(i)
    for (i = 0; i < n; ++i)
      y[incy*i] = a * x[incx*i];

  }

  template<class T, class itype=int>
  inline void 
  yeqax(const itype n, const T a, const T *x, T *y)
  {
    itype i;
   #pragma omp parallel for default(none)		\
     shared(x,y) private(i)
    for (i = 0; i < n; ++i)
      y[i] = a*x[i];
  }

  template<class T, class itype=int>
  inline void
  yeqaxpy(const itype n, const T a, const T *x, const itype incx,
	  T *y, const itype incy)
  {
    itype i;
   #pragma omp parallel for default(none)		\
     shared(x,y) private(i)
    for (i = 0; i < n; ++i)
      y[incy*i] = a * x[incx*i] * y[incy*i];
  }

  template<class T, class itype=int>
  inline void
  yeqaxpy(const itype n, const T a, const T *x, T *y)
  {
    itype i;
   #pragma omp parallel for default(none)		\
     shared(x,y) private(i)
    for (i = 0; i < n; ++i)
      y[i] = a * x[i] * y[i];
  }
  
  
  
  
  template<class T, class itype=int>
  inline void 
  xeqax(const itype n, const T a, T *x, const itype incx)
  {
    itype i;
   #pragma omp parallel for default(none)		\
     shared(x) private(i)
    for (i = 0; i < n; ++i)
      x[incx*i] *= a;
  }
  
  template<class T, class itype=int>
  inline void 
  xeqax(const itype n, const T a, T *x)
  {
    itype i;
   #pragma omp parallel for default(none)		\
     shared(x) private(i)
    for (i = 0; i < n; ++i)
      x[i] *= a;
  }
  
  template<class T, class itype=int>
  inline void 
  yeqaxy(const itype n, const T a, const T *x, const itype incx,
	 T *y, const itype incy)
  {
    itype i;
   #pragma omp parallel for default(none)		\
     shared(x,y) private(i)
    for (i = 0; i < n; ++i)
      y[incy*i] *= a * x[incx*i]; 
  }
  
  template<class T, class itype=int>
  inline void 
  yeqaxy(const itype n, const T a, const T *x, T *y)
  {
    itype i;
   #pragma omp parallel for default(none)		\
     shared(x,y) private(i)
    for (i = 0; i < n; ++i)
      y[i] *- a * x[i];
  }
  
  template<class T, class itype=int>
  inline void 
  zeqaxy(const itype n, const T a, const T* x, const itype incx, const T* y,
	 const itype incy, T* z, const itype incz)
  {
    itype i;
   #pragma omp parallel for default(none)		\
     shared(x,y,z) private(i)
    for (i = 0; i < n; ++i)
      z[incz*i] = a * x[incx*i] * y[incy*i];
  }
  
  template<class T, class itype=int>
  inline void 
  zeqaxy(const itype n, const T a, const T* x, const T* y, T* z)
  {
    itype i;
   #pragma omp parallel for default(none)		\
     shared(x,y,z) private(i)
    for (i = 0; i < n; ++i)
      z[i] = a * x[i] * y[i];
  }
  
  
  template<class T, class itype=int>
  inline void 
  zeqaxypz(const itype n, const T a, const T* x, const itype incx,
	   const T* y, const itype incy, T* z, const itype incz)
  {
    itype i;
   #pragma omp parallel for default(none)		\
     shared(x,y,z) private(i)
    for (i = 0; i < n; ++i)
      z[incz*i] += a * x[incx*i] * y[incy*i];
  }
  
  template<class T, class itype=int>
  inline void 
  zeqaxypz(const itype n, const T a, const T* x, const T* y, T* z)
  {
    itype i;
   #pragma omp parallel for default(none)		\
     shared(x,y,z) private(i)
    for (i = 0; i < n; ++i)
      z[i] += a * x[i] * y[i];
  }

  
  template<class T, class itype=int>
  inline void 
  zeqaxypbz(const itype n, const T a, const T* x, const itype incx,
	    const T* y, const itype incy, const T b, T* z, const itype incz)
  {
    itype i;
   #pragma omp parallel for default(none)		\
     shared(x,y, z) private(i)
    for (i = 0; i < n; ++i)
      z[incz*i] = a * x[incx*i] * y[incy*i] + b * z[incz*i];
  }
  


  template<class T, class itype=int>
  inline void 
  zeqaxypbz(const itype n, const T a, const T* x, const T* y, const T b, T* z)
  {
    itype i;
   #pragma omp parallel for default(none)		\
     shared(x,y,z) private(i)
    for (i = 0; i < n; ++i)
	z[i] = a * x[i] * y[i] + b * z[i];
  }

  template<class T, class itype=int>
  inline void 
  zeqawxypz(const itype n, const T a, const T* w, itype incw,
	    const T* x, const itype incx, const T* y, const itype incy,
	    T* z, const itype incz)
  {
    itype i;
   #pragma omp parallel for default(none)		\
     shared(w,x,y,z) private(i)
    for (i = 0; i < n; ++i)
      z[incz*i] += a * w[incw*i] * x[incx*i] * y[incy*i];
  }

  template<class T, class itype=int>
  inline void 
  zeqawxypz(const itype n, const T a, const T* w, const T* x,const T* y, T* z)
  {
    itype i;
   #pragma omp parallel for default(none)		\
     shared(w,x,y,z) private(i)
    for (i = 0; i < n; ++i)
      z[i] += a * w[i] * x[i] * y[i];
      }
  }
  
  template<class T, class itype=int>
  inline void
  zeqvwpxy(const itype n, const T* v, const itype incv,
	   const T* w, const itype incw,
	   const T* x, const itype incx, const T* y, const itype incy,
	   T* z, const itype incz)
  {
    itype i;
   #pragma omp parallel for default(none)		\
     shared(v,w,x,y,z) private(i)
    for (i = 0; i < n; ++i)
	z[incz*i] = v[incv*i]* w[incw*i] + x[incx*i] * y[incy*i];
  }
  
  template<class T, class itype=int>
  inline void
  zeqvwpxy(const itype n, const T* v, const T* w, const T* x, const T* y, T* z)
  {
    itype i;
   #pragma omp parallel for default(none)		\
     shared(v,w,x,y,z) private(i)
    for (i = 0; i < n; ++i)
      z[i] = v[i] * w[i] + x[i] * y[i];
    
  }
  
  
  template<class T, class itype=int>
  inline T 
  sum_ax(const itype n, const T a, const T *x, const itype incx)
  {
    T sum = 0.0;
    itype i;
   #pragma omp parallel for default(none)\
     shared(x) private(i) reduction(+:sum)
    for (i = 0; i < n; ++i)
      sum += x[incx*i];
    return a*sum;
  }

  template<class T, class itype=int>
  inline T 
  sum_ax(const itype n, const T a, const T *x)
  {
    T sum = 0.0;
    itype i;
   #pragma omp parallel for default(none)\
     shared(x) private(i) reduction(+:sum)
    for (i = 0; i < n; ++i)
      sum += x[i];
    return a*sum;
  }
  
  
  template<class T, class itype=int>
  inline T 
  sum_axy(const itype n, const T a, const T* x, const itype incx,
	  T* y, const itype incy)
  {
    T sum = 0.0;
    itype i;
   #pragma omp parallel for default(none)\
     shared(x,y) private(i) reduction(+:sum)
    for (i = 0; i < n; ++i)
      sum += x[incx*i] * y[incy*i];
    return a*sum;
  }
  
  
  template<class T, class itype=int>
  inline T 
  sum_axy(const itype n, const T a, const T* x, const T* y)
  {
    T sum = 0.0;
    itype i;
   #pragma omp parallel for default(none)\
     shared(x,y) private(i) reduction(+:sum)
    for (i = 0; i < n; ++i)
      sum += x[i] * y[i];
    return a*sum;
  }
  
  
  template<class T, class itype=int>
  inline T 
  sum_axyz(const itype n, const T a, const T* x, const itype incx,const  T* y, const itype incy, const T* z, const itype incz)
  {
    T sum = 0.0;
    itype i;
   #pragma omp parallel for default(none)\
     shared(x,y,z) private(i) reduction(+:sum)
    for (i = 0; i < n; ++i)
      sum += x[incx*i] * y[incy*i] * z[incz*i];
    return a*sum;
  }
  
  template<class T, class itype=int>
  inline T
  sum_axyz(const itype n, const T a, const T* x, const T* y, const T* z)
  {
    T sum = 0.0;
    itype i;
   #pragma omp parallel for default(none)\
     shared(x,y,z) private(i) reduction(+:sum)
    for (i = 0; i < n; ++i)
      sum += x[i] * y[i] * z[i];
    return a*sum;
  }
  
  template<class T, class itype=int>
  inline T
  sum_awxyz(const itype n, const T a, const T* w, const itype incw,
	    const T* x, const itype incx, const T* y, const itype incy,
	    const T* z, itype incz)
  {
    T sum = 0.0;
    itype i;
   #pragma omp parallel for default(none)\
     shared(x,y,z) private(i) reduction(+:sum)
    for (i = 0; i < n; ++i)
      sum += w[incw*i] * x[incx*i] * y[incy*i] * z[incz*i];
    return a*sum;
  }
  
  
  template<class T, class itype=int>
  inline T 
  sum_awxyz(const itype n, const T a, const T* w, const T* x, const T* y, const T* z)
  {
    T sum = 0.0;
    itype i;
   #pragma omp parallel for default(none)\
     shared(w,x,y,z) private(i) reduction(+:sum)
    for (i = 0; i < n; ++i)
      sum += w[i] * x[i] * y[i] * z[i];
    return a*sum;
  }
  
  
  
  template<class T, class itype=int>
  inline void 
  BeqDxA(const itype n, const itype m, const T *x, const itype incx,
      const T *A, const itype lda, T *B, const itype ldb)
  {

    itype i;
   #pragma omp parallel for default(none)\
     shared(x,A,B) private(i)
    for (i = 0; i < n; ++i)
      yeqax(m, x[incx*i], A + i*lda, 1, B+ldb*i, 1);
  }
  
  template<class T, class itype=int>
  inline void 
  BeqADx(const itype n, const itype m, const T *A, const itype lda,
      const T *x, const itype incx, T *B, const itype ldb)
  {
    itype i;
   #pragma omp parallel for default(none)\
     shared(x,A,B) private(i)
    for (i = 0; i < n; ++i)
      zeqaxy(m, 1.0, x, incx, A + lda*i, 1, B + ldb*i, 1);
  }
  
  template<class T, class itype=int>
  inline void 
  BeqDxADy(const itype n, const itype m, T *x, const itype incx,
      const T *A, const itype lda, const T *y, const itype incy,
      T *B, const itype ldb)
  {
    itype i;
   #pragma omp parallel for default(none)\
     shared(x,A,B) private(i)
    for (i = 0; i < n; ++i)
      zeqaxy(m, x[incx*i], y, incy, A + lda*i, 1, B + ldb*i, 1);
  }
  
  template<class T, class itype=int>
  inline
  void y_eq_trans_x(const itype n , const T *x, itype *trans, T *y){

    itype i;
    
#pragma omp parallel for default(none)\
     shared(x,y,trans) private(i)
    for (i = 0; i < n; ++i)
      y[i] = x[trans[i]];
  }
  
  template<class T, class itype=int>
  inline
  void yeqmapx(const itype n, const T *x, const itype *map, T *y)
  {
    itype i;
    
#pragma omp parallel for default(none)\
     shared(x,y,map) private(i)
    for (i = 0; i < n; ++i)
      y[i] = x[map[i]];
  }

  template<class T, class itype=int>
  inline
  void yeqmapx(const itype n, const T *x, const itype *map, T *y, itype incy)
  {
    itype i;
    
#pragma omp parallel for default(none)\
     shared(x,y,map) private(i)
    for (i = 0; i < n; ++i)
      y[incy*i] = x[map[i]];
	y += incy;
  }

  template<class T, class itype=int>
  inline
  void mapyeqx(const itype n, T *x, itype *rmap, T *y)
  {
    while (n--)
      {
	y[*rmap++] = *x++;
      }
  }
  

  template<class T, class itype=int>
  inline
  void mapyeqx(const itype n, const T *x, const itype incx,
	       const itype *rmap, T *y)
  {
    itype i;
    
#pragma omp parallel for default(none)\
     shared(x,y,rmap) private(i)
    for (i = 0; i < n; ++i)
      y[rmap[i]] = x[incx*i];
 
}


#endif // _MMSL_VEC__HPP_
